Skip to content

feat: add SemilinearODEFunction and SemilinearODEProblem #3739

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom
Open
wants to merge 10 commits into
base: master
Choose a base branch
from
Open

feat: add SemilinearODEFunction and SemilinearODEProblem #3739

wants to merge 10 commits into from

Conversation

AayushSabharwal
Copy link
Member

Generates

f1 = B * x2 + C
f2 = A * x

Where x2 is the upper triangle of x * x' as a linear vector with unnecessary elements zeroed. This is what Symbolics.semiquadratic_form returns.

In terms of codegen, the IIP form uses a DiffCache to avoid allocating a new buffer for x2 each time (and to allow passing duals) and generates the equivalent of

x2 = view(PreallocationTools.get_tmp(...), ...)
x2[1] = ...
x2[2] = ...
...
out[1] = C[1]
out[2] = C[2]
...
mul!(C, B, x2, true, true)

The jacobian for f1 is B * x2' + C' where x2' is the jacobian of x2 and C' that of C. This again uses the DiffCache and generates code with the same structure as above. For sparse forms, it is basically identical except the array from the DiffCache is used as the nzvals of a SparseMatrixCSC.

Still needs a little more implementation for the sparse form, tests and docs.

@AayushSabharwal AayushSabharwal marked this pull request as draft June 16, 2025 08:18
@AayushSabharwal
Copy link
Member Author

No point running CI right now if no tests or docs are changed.

@AayushSabharwal AayushSabharwal force-pushed the as/semilinear-odeprob branch from b9dc547 to bbda8b1 Compare June 25, 2025 14:47
@AayushSabharwal AayushSabharwal marked this pull request as ready for review June 25, 2025 14:47
@AayushSabharwal AayushSabharwal force-pushed the as/semilinear-odeprob branch from bbda8b1 to 645a720 Compare July 9, 2025 06:24
@AayushSabharwal AayushSabharwal force-pushed the as/semilinear-odeprob branch from 645a720 to 6b27ebc Compare July 16, 2025 08:27
@ctessum
Copy link
Contributor

ctessum commented Aug 11, 2025

@AayushSabharwal thanks for this!

I tried it with the following example (from here):

using Pkg
Pkg.add(url="https://github.com/AayushSabharwal/ModelingToolkit.jl", rev="as/semilinear-odeprob")

using ModelingToolkit, MethodOfLines
using OrdinaryDiffEqSDIRK
t = ModelingToolkit.t_nounits
D = ModelingToolkit.D_nounits
using DomainSets
using MethodOfLines
using Plots


x_min = y_min = t_min = 0.0
x_max = y_max = 1.0
t_max = 10.0

N = 10

dx = (x_max-x_min)/N
dy = (y_max-y_min)/N

islocation(x, y) = x > x_max / 2 - dx && x < x_max / 2 + dx &&
    y > y_max / 10 - dx &&
    y < y_max / 10 + dx
emission(x, y) = ifelse(islocation(x, y), 10, 0)
@register_symbolic emission(x, y)

θ(x,y) = atan(y.-0.5, x.-0.5)
u(x,y) = -sin(θ(x,y))
v(x,y) = cos(θ(x,y))

@parameters T=293.15 x y

@constants k=0.1 T_0=300.0
@variables A(..) B(..)
Dx = Differential(x)
Dy = Differential(y)
advect(var) = -u(x,y)*Dx(var) - v(x,y)*Dy(var)
eqs = [
    D(A(x,y,t)) ~ advect(A(x, y, t)) + emission(x, y) - k*exp(T/T_0)*A(x,y,t),
    D(B(x,y,t)) ~ advect(B(x, y, t)) + k*exp(T/T_0)*A(x,y,t),
]


domain = [
    x ∈ Interval(x_min, x_max),
    y ∈ Interval(y_min, y_max),
    t ∈ Interval(t_min, t_max),
]

bcs = [A(x,y,t_min) ~ 0.0,
       A(x_min,y,t) ~ A(x_max,y,t),
       A(x,y_min,t) ~ A(x,y_max,t),

       B(x,y,t_min) ~ 0.0,
       B(x_min,y,t) ~ B(x_max,y,t),
       B(x,y_min,t) ~ B(x,y_max,t),
]

@named pdesys = PDESystem(eqs, bcs, domain, [x,y,t], [A(x,y,t), B(x,y,t)], [T, k, T_0])

discretization = MOLFiniteDifference([x=>dx, y=>dy], t);

odesys, tspan = symbolic_discretize(pdesys, discretization)

odesys_complete = complete(odesys)

# Gives error: ERROR: Encountered operator `Differential(t)` which has different independent variable than the one used in the system `nothing`.
prob = SemilinearODEProblem(odesys_complete, [], tspan; stiff_linear = false,
    stiff_quadratic = true, stiff_C = true)

prob = ODEProblem(odesys_complete, [], tspan) # This gives the same error.

Here is the stack trace:

ERROR: Encountered operator `Differential(t)` which has different independent variable than the one used in the system `nothing`.

Context:
0.7071067811865476(-10.0(A(t))[11, 2] + 10.0(A(t))[2, 2]) - 0.7071067811865475(-10.0(A(t))[2, 2] + 10.0(A(t))[2, 3]) + Differential(t)((A(t))[2, 2]) + (A(t))[2, 2]*k*exp(T / T_0)
Stacktrace:
  [1] validate_operator(op::Differential, args::SymbolicUtils.SmallVec{…}, iv::Nothing; context::SymbolicUtils.BasicSymbolic{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/QdAET/src/utils.jl:541
  [2] collect_vars!(unknowns::OrderedCollections.OrderedSet{…}, parameters::OrderedCollections.OrderedSet{…}, expr::SymbolicUtils.BasicSymbolic{…}, iv::Nothing; depth::Int64, op::Type)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/QdAET/src/utils.jl:599
  [3] collect_vars!(unknowns::OrderedCollections.OrderedSet{…}, parameters::OrderedCollections.OrderedSet{…}, eq::Equation, iv::Nothing; depth::Int64, op::Type)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/QdAET/src/utils.jl:620
  [4] collect_scoped_vars!(unknowns::OrderedCollections.OrderedSet{…}, parameters::OrderedCollections.OrderedSet{…}, sys::System, iv::Nothing; depth::Int64, op::Type)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/QdAET/src/utils.jl:507
  [5] collect_scoped_vars!
    @ ~/.julia/packages/ModelingToolkit/QdAET/src/utils.jl:500 [inlined]
  [6] discover_globalscoped(sys::System)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/QdAET/src/systems/abstractsystem.jl:608
  [7] complete(sys::System; split::Bool, flatten::Bool, add_initial_parameters::Bool)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/QdAET/src/systems/abstractsystem.jl:631
  [8] complete
    @ ~/.julia/packages/ModelingToolkit/QdAET/src/systems/abstractsystem.jl:629 [inlined]
  [9] macro expansion
    @ ~/.julia/packages/Setfield/ZezIj/src/sugar.jl:198 [inlined]
 [10] mtkcompile(sys::System; additional_passes::Vector{…}, simplify::Bool, split::Bool, allow_symbolic::Bool, allow_parameter::Bool, conservative::Bool, fully_determined::Nothing, inputs::Vector{…}, outputs::Vector{…}, disturbance_inputs::Vector{…}, kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/QdAET/src/systems/systems.jl:51
 [11] mtkcompile
    @ ~/.julia/packages/ModelingToolkit/QdAET/src/systems/systems.jl:30 [inlined]
 [12] InitializationProblem{…}(sys::System, t::Float64, op::Dict{…}; guesses::Dict{…}, check_length::Bool, warn_initialize_determined::Bool, initialization_eqs::Vector{…}, fully_determined::Nothing, check_units::Bool, use_scc::Bool, allow_incomplete::Bool, algebraic_only::Bool, time_dependent_init::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/QdAET/src/problems/initializationproblem.jl:66
 [13] InitializationProblem
    @ ~/.julia/packages/ModelingToolkit/QdAET/src/problems/initializationproblem.jl:20 [inlined]
 [14] #_#1116
    @ ./none:0 [inlined]
 [15] maybe_build_initialization_problem(sys::System, iip::Bool, op::Dict{…}, t::Float64, defs::Dict{…}, guesses::Dict{…}, missing_unknowns::Set{…}; implicit_dae::Bool, time_dependent_init::Bool, u0_constructor::Function, p_constructor::Function, floatT::Type, initialization_eqs::Vector{…}, use_scc::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/QdAET/src/systems/problem_utils.jl:1125
 [16] process_SciMLProblem(constructor::Type, sys::System, op::Dict{…}; build_initializeprob::Bool, implicit_dae::Bool, t::Float64, guesses::Dict{…}, warn_initialize_determined::Bool, initialization_eqs::Vector{…}, eval_expression::Bool, eval_module::Module, fully_determined::Nothing, check_initialization_units::Bool, u0_eltype::Nothing, tofloat::Bool, u0_constructor::typeof(identity), p_constructor::typeof(identity), check_length::Bool, symbolic_u0::Bool, warn_cyclic_dependency::Bool, circular_dependency_max_cycle_length::Int64, circular_dependency_max_cycles::Int64, substitution_limit::Int64, use_scc::Bool, time_dependent_init::Bool, algebraic_only::Bool, allow_incomplete::Bool, is_initializeprob::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/QdAET/src/systems/problem_utils.jl:1384
 [17] SemilinearODEProblem{…}(sys::System, op::Vector{…}, tspan::Tuple{…}; check_compatibility::Bool, u0_eltype::Nothing, expression::Type, callback::Nothing, sparse::Bool, stiff_linear::Bool, stiff_quadratic::Bool, stiff_nonlinear::Bool, jac::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/QdAET/src/problems/odeproblem.jl:209
 [18] SemilinearODEProblem
    @ ~/.julia/packages/ModelingToolkit/QdAET/src/problems/odeproblem.jl:162 [inlined]
 [19] SemilinearODEProblem
    @ ./none:0 [inlined]
 [20] #SemilinearODEProblem#1020
    @ ./none:0 [inlined]

@ctessum
Copy link
Contributor

ctessum commented Aug 11, 2025

Update, if I do this:

odesys_complete = mtkcompile(odesys)

then it no longer gives that error. Now I get the error:

ERROR: ArgumentError: All of `A`, `B` and `C` cannot be non-stiff at the same time.

I guess this is because, in this example, both the chemistry terms and the advection terms are linear?
In general what we want here is for the chemistry terms to be considered stiff and the advection terms to be considered non-stiff.

Would it make sense for there to be an option for the user to supply functions like linear_stiff(term)::Bool, quadratic_stiff(term)::Bool, etc, to give the user control of which terms get put into which functions?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers
No reviews
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.
2 participants
@AayushSabharwal @ctessum